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We review the local Monte Carlo dynamics and Swendsen-Wang cluster algo- 
rithm. We introduce and analyze a new Monte Carlo dynamics known as tran- 
sitional Monte Carlo. The transitional Monte Carlo algorithm samples energy 
probability distribution P(E) with a transition matrix obtained from single-spin- 
flip dynamics. We analyze the relaxation dynamics master equation, 



dP{E,t) 
It 



^2T(E,E')P(E',t), 



associated with fsing model in d dimensions. In one dimension, we obtain an 
exact solution. We show in all dimensions in the continuum limit the dynamics 
is governed by the partial differential equation 

dP _ cfP J)P_ 

dt' dx 2 dx 

where x and t' are rescaled energy deviation from the equilibrium value and 
rescaled time, respectively. This equation is readily solved. Thus, we have a 
complete understanding of the dynamics. 



I. MONTE CARLO METHOD 

Monte Carlo method in the most basic application is to perform numerical integration 
of very high dimensional integrals or to compute averages of a given probability distribution. 
In this respect, the method generates a sequence of states X , Xi, X 2l ■ ■ . , by a transition 
probability T(X — > X') = P(X'\X). This is also known as Markov chain Monte Carlo in the 
statistics community. The state X can be the set of all the coordinates of the particles in a 
fluid systems, or the values of spins at lattice points for a classical spin system. The new state 
X' is generated according to probability P(X'\X) given that the previous state is X. P{X'\X) 
is usually a simple distribution which can be sampled directly. 

We have a lot of freedom in choosing the transition probability T(X — > X'). If we want that 
the distribution of X of the generated sequence follows P eq (X), it is sufficient to require that 

P eq (X) T(X -» X 1 ) = P eq (X') T(X' -» X). (1) 

This is known as detailed balance condition. Subject to some general constraints of ergodicity 
which can be satisfied easily, the state X has the unique equilibrium distribution in the large 
step limit. 

Well-known application of Monte Carlo method is to generate Boltzmann distribution 
e -H/k B T J % j n statistical mechanics by Metropolis algorithm, where the transition probabil- 
ity is chosen as 
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T(X^X') = W(X^X')min(l ) ^^j ) X + X\ 



(2) 



where W(X — » X') = W(X' — » X) is the probability of proposing X' as the next state given 
that the current state is X. The next factor min(- ■ ■) is the probability that such a move is 
accepted. 

Averages of thermodynamic variables are computed according to 

1 M 

(A) = £ A(X)P eq (X) » - £ (3) 

The weighted sum (or integral) is replaced by an arithmetic average. 

Monte Carlo method is intrinsically an approximate method. Error in Monte Carlo evaluation 
can be estimated from 

M«-^=, (4) 

where a \ is the variance of distribution of A, and r is called correlation time. The value 
t > 1 is more precisely defined by so-called integrated correlation time. But in many cases it 
is equivalent to the exponential correlation time defined by the equation: 

(A(t )A(t + t)) - (A(t )) 2 oc e- l l\ (5) 

The exponential correlation time r also characterizes the speed at which arbitrary initial prob- 
ability distribution P(X) converges to P eq (X): 

P{X, t) » P eq (X) + Pi{X)e-V T + ■■■. (6) 



II. SOME MONTE CARLO DYNAMICS 



The above theory is quite general. We give some more concrete examples and points out 
some interesting Monte Carlo dynamics. 



A. Single-spin-flip Glauber dynamics 



We consider a simple classical spin model, the Ising model, as an example of local Monte 
Carlo dynamics. The model is defined by the energy function 

H(a) = -jJ2wji ff< = ±l. (7) 

(i,3) 



The spins take only two possible values and live on the sites of a lattice, for example on a 
square lattice. The total energy is a sum of interactions between nearest-neighbor sites. A 
Monte Carlo move consists of picking a site at random, and flipping the spin with probability 
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where ctq is the spin value before the flip, and <ii is the sum of spins of the nearest neighbors. 
Another popular choice is the Metropolis rate min(l, exp(— AH/ksT)) where AH is the energy 
increase due to flip. In a computer implementation, the value w is compared with a uniformly 
distributed random number r between and 1. If r < w the flip is performed; otherwise, the 
old configuration is retained and is also counted as one Monte Carlo move. One Monte Carlo 
step is a unit (Monte Carlo) time, and is usually defined as TV Monte Carlo moves, where N is 
the number of spins of the system. 

The local Monte Carlo dynamics has some common features: (1) the algorithm is extremely 
general. It can be applied to any classical model. (2) Each move involves 0(1) operation and 
O(l) degrees of freedom. (3) At the second-order phase transition critical point, the dynamics 
becomes very slow. This is known as dynamical critical slowing down, characterized by the fact 
that 

rcx<f, 2»2. (9) 

where £ oc \T — T c \~ v is called correlation length which will diverge at the critical point. On 
finite system, £ is replaced by the system linear size L. The value z is the dynamical critical 
exponent. Its value is slightly greater than 2 for a large class of models with order-parameter 
nonconserving dynamics, such as the single-spin-flip dynamics discussed above Q. 

This last feature hampers the effective use of local Monte Carlo algorithms. It is the nonlocal 
algorithms that come to rescue. 



B. Nonlocal dynamics — cluster algorithms 

Swendsen-Wang algorithm |||(| is one of the first nonlocal Monte Carlo algorithms which have 
very different dynamical characteristics. The algorithm uses a mapping from Ising model to a 
type of percolation model. Each Monte Carlo step consists of putting a bond with probability 

p{ai,<Tj) = 1 - exp(- JipiOj + l)/k B T) (10) 

between each pair of the nearest neighbors, by ignoring the spins and looking only at the bonds, 
we obtain a percolation configurations of bonds ]8| . A new spin configuration is obtained by 
assigning to each cluster, including isolated sites, a random sign +1 or —1 with equal probability. 

In the Swendsen-Wang algorithm, we generated many clusters and then flipped these clusters. 
Wolff algorithm |?J is a variation on the way clusters are flipped. One picks a site at random, 
and then generates one single cluster by growing a cluster from the seed. The neighbors of the 
seed site will also belong to the cluster if the spins are in parallel and a random number is less 
than p=l — e - 2J / kBT , xhat is, the neighboring site will be in the same cluster as the seed site 
with probability p if the spins have the same sign. If the spins are different, neighboring site will 
never belong to the cluster. Neighbors of each new site in the cluster are tested for membership. 
This testing of membership is performed on pair of sites (forming a nearest neighbor bond) not 
more than once. The recursive process will eventually terminate. The spins in the cluster are 
turned over with probability 1. 

The following is a fairly elegant way of implementing the Wolff algorithm in C. The function 
is the core part which performs a Wolff single cluster flip. This function is recursive. The array 
for spins s [ ] , percolation probability P, and coordination number Z (the number of neighbors) 
are passed globally. The first argument i of the flip function is the site to be flipped, the 
second argument sO is the spin of the cluster before flipping. The function neighbor returns 
an array of neighbors of the current site. The function drand48() returns a random number 
unformly distributed between and 1. 



3 



void flip(int i, int sO) 
{ 

int j , nn [Z] ; 

s [i] = - sO; flip the spin immediately 

neighbor (i, nn) ; find nearest neighbor of i 

for(j =0; j < Z; ++j) flip the neighbor if 

if(s0 == s[nn[j]] && drand48() < P) spins are equal and 

flip(nn[j] , sO) ; random number is smaller than p 

} 

Some of the salient features of cluster algorithms: (1) algorithm is applicable to models 
containing Ising symmetry. (2) Computational complexity is still of O(l) per spin per Monte 
Carlo step. (3) Much reduced critical slowing down, i.e., r cx £ Zs ™. The dynamical critical 
exponent z sw is 0, 0.3, 0.5, 1, in dimensions 1, 2, 3, and greater or equal to 4, respectively. In 
addition, Li and Sokal Q showed that r > ac for some constant a, and c is the specific heat. 
For some of the lastest developments on nonlocal algorithms, see 10 - 13| . 



III. TRANSITIONAL MONTE CARLO DYNAMICS 



The transitional Monte Carlo dynamics [|14j is a new dynamics with the following interesting 
features: (1) It is a constrained random walk in energy space. (2) The transitional rates are 
derived from single-spin-flip dynamics. (3) It has a fast dynamics, r cx c, and (4) it suggests a 
different histogram reweighting method. 



A. What is transitional Monte Carlo 



The transitional Monte Carlo is related to the single-spin-flip dynamics in the following sense 
but it has a totally different dynamics. A single-spin- flip Glauber dynamics of the Ising model 
is described by 



{*'} 



N 

^ [-Wi(cr») + Wi(-(Ji)Fi P{cr,t), 



(11) 



where N is the total number of spins, and w is given by Eq. (|J), and Fi is a flip operator such 
that FjP(. . . , <7j, . . .) = P(. . . , —(Ti, . . .). Transitional Monte Carlo dynamics is defined by 



^§^-=J2T(E,E')P(E', t ), 



(12) 



(13) 



where P(E, t) is the probability of having energy E at time t, and 

t ^e') = ^- £ £ ivx), 

' H(a)=E H{a')=E' 

where n(E) is the degeneracy of the states. We can not derive Eq. ( jl2) from Eq. ( pd| ) 
in general, the "derivation" is valid only at equilibrium when P{E) = z2,H{a)=E^'i tJ ) = 
n(E)exp(~E/k B T). 
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The transition matrix T(E, E') has some general properties: (1) The matrix is banded alone 
the diagonal. (2) The column sum is zero, ^ E T{E,E') — 0, due to the conservation of total 
probability. (3) ^2 E ,T{E 1 E')P eq {E') — 0, due to existence of equilibrium distribution. (4) 
The transition rate satisfies detailed balance condition, T(E' ', E)P eq (E) = T(E, E')P eq (E'). 
The eigenvalues of T(E,E') are real and Aft < 0. The significance of the eigenspectrum is 
that the general solution of the master equation can be written as P(E, t) — J2k a k(E) exp(Afct). 



B. Computer realization of transitional Monte Carlo and reweighting techniques 

The transitional Monte Carlo dynamics can be implemented on computer in at least two 
different ways, we'll call them algorithm A and B. 

Algorithm A 

1. Do sufficient number of constant energy (microcanonical) Monte Carlo steps, so that 
the final configuration is totally uncorrelated with the initial configuration. This step is 
equivalent to pick a state a at random from all states with energy E. 

2. Do one canonical Monte Carlo move by picking a site at random. 

Clearly, this algorithm is not very efficient computationally, due to step 1. However, it will be 
helpful in understanding the dynamics. 

Algorithm B 

• A direct implementation of the master Eq. (p"2|), i.e., a random walk in energy with a 
transition rate T(E, E'). 

Other possibility is to solve the equation on computer. Then in this method and algorithm 
B, we need to know T(E, E') explicitly, this can be done numerically by Monte Carlo sampling, 
from 

T(E + AE, E) = w(AE)(N(a, AE)) H{a)=E (14) 

and w{AE) = ±(l - tanh(A£/(2fc B T)). N(a,AE) is the number of cases that energy is 
changed by AE from E for the N possible single-spin flips. 

Note that computation of (N(cr, AE)} h , ._ e can be done with any sampling technique which 
ensures equal probability for equal energy. It is a kind of "combinatorial" number independent 
of the spin flip rates and in particular, independent of the temperature. Thus the transition 
matrix can be formed with any temperature. The equilibrium distribution and thus the density 
of states n(E) = P eq {E) exp(E/ksT) is obtained by solving 

Y J T{E,E')P eq {E') = Q. (15) 

E' 

The above scheme is similar in spirit to the histogram method of Ferrenberg and Swendsen |l5f| , 
and the method has a close connection with, but different from the broad histogram of Olivcira 
et al (H). 
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C. Exact results in transitional Monte Carlo dynamics 



We have more or less a complete understanding of the transitional Monte Carlo dynamics 
through exact results in limiting cases. The transition matrix T(E, E') can be computed exactly 
in one-dimensional chain of length L (with periodic boundary condition) , by some combinatorial 
consideration, as 

T k , k+ i = i j^j '-{I + 7), (16) 

(L-2k)(L-2k-l) 
Tk+i,k= _ ^ 1 1 -!), i 1 ') 

where 7 = tanh(2 J/ksT). The diagonal term is computed from the relation 

T k -x, k + T kik + T k+1>k = 0, (18) 

and the rest of the elements T k<k i = Q if \k — k'\ > 1. The integer k = 0, 1, 2, ... , [L/2\ is 
related to energy by E = —LJ + 4k. While the eigen spectrum at temperature T — can be 
computed exactly as = —2(k + l)(2k + l)/(£ — 1), the eigenvalues at T > is obtained only 
numerically. The most important feature is that r oc L, given an unusual dynamical critical 
exponent of z = 1. 

The dynamics in any dimensions in the large size limit obeys a linear Fokker-Planck equation: 

where t' and x are properly scaled time and energy. 

E-u N 



{Nc') 1 / 2 ' 



u Q N = E, (20) 



and t' = bt with 



b= Urn — !— Yt(E,E')(E' - E) 2 , (21) 
n^oo 2c'N ^— ' y ' v y 

E' 

where uo is the average energy per spin and d = ksT 2 c is the reduced specific heat per spin. 
The major consequence of this result is that the relaxation times are r n — ad/n, n = 1, 2, 3, • • •, 
with some constant a. 

We can cast the equation in the form of a continuity equation, 
dP dj dP 

_ + -2=0, with J— 5 -«P. (22) 

There are two competing effects in the current; —dP/dx is the usual diffusion, while — xP is a 
counter drift, j = produces the equilibrium distribution P eq (x) oc exp(— x 2 /2). 

With a change of variable P(x,t') = exp(— x 2 /4)<&{x, t'), the Eq. ( ^9|) is transformed into a 
one-dimensional quantum harmonic oscillator equation, with a general solution of the form 

00 

J2 c neM-nbt-x 2 /2)H n (x/V2) 7 (23) 

71 = 



G 



where H n is Hermite polynomials. 

We sketch the derivation of the continuum limit equation, which is known as Q expansion IlTj] . 
Starting from the master equation, Eq. (|l2|), we introduce the new variable x which describes 
the scaled deviations from equilibrium, 

E = E + 5E = u Q N + x(Nc') 1/2 . (24) 

The function P(E, t) is written in terms of x, and P(E', t) — > P(x + Sx, t) = P(x, t) + ||&r + 

^^S-Sx 2 + • • •. For the matrix T(E,E'), we assume that the changes along the diagonal are 
smooth and can be expanded, but across diagonals the changes are still discrete. For T(E, E') 
we also expand around x — 0; we can show that such an expansion is also an expansion in 
power of N~ x / 2 . Expanding all the relevant terms in powers of A -1 / 2 , the leading terms of 
order N° give the desired equation. The rest of the correction terms go to zero in the large size 
limit N — > oo. We used some of the general properties of T(E, E') to simplify the equation. 



D. Simple pictures of the dynamics 



The exact results can be interpreted with intuitive pictures. First, we consider the result of 
r cx L in one dimension as T — > 0. At sufficiently low temperatures with a correlation length £ 
comparable with the system size L, only the ground state (all spins up or down) and the first 
excited states (with a kink pair) are important. Let's consider the time scale for Eq — > E\, a 
spin with opposite sign is created with probability exp(— AK) from Boltzmann weight, where 
K = J/(fcsT), in each of the canonical move. Thus 

exp(4AQ e 

t cx cx — cx L. (25) 

L L 

where K is chosen such that there is about one kink pair, so that £ ~ exp(2A) ~ L. 

The reverse process, Ei — > E has also the same time scale. In this case, it is advantageous 
to use the equivalence of Algorithm A and B. Given that the system is in the state of kink pair 
(a string of up spins followed by a string of down spins, with periodic boundary condition), the 
first step of Algorithm A randomizes the locations of the kinks. The probability that two kinks 
are the nearest neighbors is X/L; the probability that this pair is chosen and destroyed by a 
spin flip is 1/L. Thus, the transitional Monte Carlo moves needed to destroy a kink pair are 
i/((l/Z)(l/L)) = L 2 . The time in terms of Monte Carlo step is then r cx L 2 / L = L. 

Similarly, the result of r cx c can be obtained by the following argument. The transitional 
Monte Carlo is a random walk constrained in the range 5E, due to the gaussian distribution 
nature of the equilibrium distribution P eq (E). The width of this distribution is related to the 
specific heat by SE 2 = cNksT 2 . Each walk changes E by 0(1). To change E by SE, we need 
5E 2 moves, invoking the theory on random walks. In units of transitional Monte Carlo steps, 

SE 2 , . 

r w a— — cx c. (26) 
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